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Abstract—Inferring causal relationships is key to data 
science. Learning causal structures in the form of directed 
acyclic graphs (DAGs) has been widely adopted for uncovering 
causal relationships, nonetheless, it is a challenging task owing 
to its exponential search space. A recent approach formulates 
the structure learning problem as a continuous constrained 
optimization task that aims to learn causal relation matrix. 
Following it are various variants, such as nonlinear variants 
which uncover nonlinear causal relationships. However, the 
nonlinear variant which considers an ordinary sparsity 
constraint as part of its optimization objective may not 
effectively eliminate spurious edges, i.e., wrongly predicted 
edges. In this paper, we investigate the defect that the ordinary 
sparsity constraint inherently contains, that is, when the model 
predicts spurious edges, the sparsity constraint may be unable 
to make the relation matrix sparse to help eliminate such 
spurious edges, and based on the theoretical and empirical 
analysis of the defect, we propose the proportional sparse 
penalty constraint PSP which replaces the ordinary sparsity 
constraint with a normalized first-order matrix norm. We then 
compare our proposed model GAE-PSP with three models to 
show considerable performance improvement. We further 
conduct three comparative experiments and show that the PSP 
can make the relation matrix sparse to eliminate spurious edges. 


Keywords—causal structure learning, proportional sparse 
penalty, continuous optimization, nonlinearity 


I. INTRODUCTION 


Discovering causal relationships among things is one of 
the core problems of science [1]. Although randomized 
controlled trial (RCT) is a classical and effective way to 
discover causal relationships, 1t may be unethical, expensive, 
or even infeasible due to experiment conditions [2]. Learning 
causal relationships from purely observational data 
circumvents these restrictions and has been flourishing in 
recent years [3]. The most common methods learn the directed 
acyclic graphs (DAGs) as causal relationships and have been 
applied to various domains such as biology [4], economics [5], 
and earth system sciences [6]. 


Causal structure learning is defined as unsupervised 
learning that given observational variables, discovers causal 
relationships among them in the form of DAG, completed 
partially directed acyclic graph (CPDAG), and others. There 
are three major classes of models, namely constraint-based, 
score-based, and those based on functional causal models 
(FCMs) [7]. Constraint-based methods, including PC [8] and 
IC [9], under faithfulness assumption, use conditional 
independent tests to recovery the Markov equivalence class 
and orient other edges. Score-based methods, such as GES 
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[10], replace conditional independent tests with the goodness 
of fit tests and try to optimize a pre-defined score function to 
find the best structure under locally consistent and 
decomposable assumptions. Algorithms based on FCMs, like 
LINGAM [11] and ANM [12], assume data are generated 
based on a structural equation model (SEM) and mutually 
independent non-Gaussian noises. However, all the methods 
have to combine the acyclicity identification step to estimate 
DAGs, thus are inefficient in practice. 


Recently, Zheng et al. [13] have formulated the structure 
learning problem as a purely continuous optimization problem 
for avoiding verifying the aforementioned combinatorial 
acyclicity constraint. In particular, they convert the acyclicity 
identification step to a part of continuous optimization 
objective and propose a score-based method which is called 
NOTEARS. The method predicts variables by variables 
themselves via the SEM and optimizes with the objective 
consists of a reconstruction error, a sparsity constraint, under 
the aforementioned acyclicity constraint. The insight of the 
method is to substitute local convergence with global 
convergence, along with SEM, to transform the pipeline 
model into the end-to-end model, which can efficiently solve 
DAG structure learning problems. Inspired by the method, 
subsequent works extend the NOTEARS framework to 
various applications [14,22], as well as nonlinear variants such 
as GAE [23], DAG-GNN [24], and NOTEARS-MLP [25]. 


However, we notice that nonlinear methods under the 
NOTEARS framework suffer from the defect: when the 
reconstruction error traps into a local optimum where the 
variable uses other non-causal variables to fit itself and thus 
introduces spurious edges, the ordinary sparsity constraint 
may not effectively help to eliminate spurious edges. Fig. 1 is 
a quick illustration of the problem: the ordinary sparsity 
constraint may not help eliminate the vertical noise line on the 
9th column of the estimated graph. More details can be seen in 
Section III. 


In this work, we investigate why the ordinary sparsity 
constraint does not work effectively. Owing to the property 
that both reconstruction error and sparsity constraint constitute 
the optimization objective, as well as the nonlinear 
transformers, whenever the reconstruction error traps into a 
local optimum, which means it adds a large number of non- 
causal variables for fitting and forms spurious edges, 
unsuitable sparsity constraint does not effectively help the 
reconstruction error escape the local optimum. Specifically, 
the optimization objective will minimize the sparsity 
constraint by amplifying the nonlinear transformers and 
decreasing the absolute maximum value of the estimated 
matrix rather than making the matrix sparse, thus the noise 


True Graph 


Estimated Graph 


Fig. 1. Visual comparison of the true relation matrix (left) and the estimated 
relation matrix (right). The matrix shows pointing edges of the 50 nodes, 
which represents the causal relationships of them. A certain point on the i" 
row and j" column means there is an edge from node i to node j. We see 
that a vertical noise line is on the 9 column of the estimated matrix, where 
the light-colored points denote the false positive predictions (spurious edges) 
except that the dark-colored one on the 23" row is true (see the red circle 
mark). The vertical noise lines often appear in experiments, and are what we 
want to mitigate. Note that prediction of an edge only relate to the lightness 
of the point, regardless of its hue, i.e., no matter red or blue. 


values of spurious edges still have a relatively large ratio to 
values of true edges. We analyze it theoretically by studying 
the optimization objective, as well as empirically prove its 
existence by conducting experiments. Inspired by that, we 
develop a constraint called proportional sparse penalty (PSP) 
to replace ordinary sparsity constraint. Our contributions are 
three-fold: 


e We theoretically and empirically explore the existence 
and the reason of the defect of ordinary sparsity 
constraint that arises from unsuitable sparsity 
constraint, which allows us to devise an effective 
approach to eliminate spurious edges with pertinence. 


e We develop a constraint called proportional sparse 
penalty (PSP) for mitigating the defect in nonlinear 
NOTEARS models, which replaces the ordinary 
sparsity constraint with a normalized first-order matrix 
norm. We further incorporate it into the GAE to 
propose a new model named GAE-PSP. 


e We demonstrate the superiority of the proposed model 
by comparing with three models and validate the 
effectiveness of the proposed constraint on three 
comparative experiments with ordinary constraint. 


II. BACKGROUND 


In this section, we first formalize the basic causal structure 
learning problem, then introduce score-based learning and 
SEM, and methods of NOTEARS class at last. 


A. Causal Structure Learning 


The causal structure learning problem is formulated as 
follows: Let X = (X,,°",Xm)! E€ R™** as a d-dimensional 
vector with m variables, and let G = (V, E) on m nodes as a 
DAG where V and EF are node set and edge set respectively. 
The DAG G represents the data generation mechanism of X, 
and We expect to learn G when given only X [24]. 


B. Score-based Learning and SEM 

Score-based algorithms learn causal structures by 
maximum the well-designed score criterion S (X i G) where 
higher score means more accurate causal structure, which is 


gig (XG), (1) 


s.t. G(W) € DAGs. 


When modeling X via linear SEM, we define W = 
[w| [Wm] E R””™ as a weighted matrix of G, where 
wi; + 0 means there is an edge from node i to node j. We 
describe the data generation mechanism of X via a linear SEM 
by X = wj X + Zz; , where Z = (Z, Zm)" E€ R™* is a 
random noise vector of any distribution types. We can rewrite 
the linear SEM as 


X=W'X+4+Z. (2) 


SEM is one of popular models used in applied statistics 
which aims to uncover inherent structural relationships 
between variables [26]. In this paper, we only use its structural 
model definition: 


n = Bn +I% +, (3) 


where 7 and € represent endogenous and exogenous variables 
respectively, ¢ is the error term, B and I are coefficient 
matrices. Note that exogenous variables affect endogenous 
variables but the reverse is not true. Common statistical 
methods for estimating SEM are unweighted least squares 
(ULS), generalized least squares (GLS) and maximum 
likelihood estimation (MLE). 


C. The NOTEARS Method 


Zheng et al. [13] are generally considered as the first to 
transform the combinatorial score-based structure search 
problem into a continuous optimization problem. They focus 
on linear SEM and the mean squared error (MSE) aka 
reconstruction error which is 


L(W;X ~ + Wy wx? 4 
( J ) F I (4) 


and add sparsity constraint like £} -regularization ||W ||} to 
learn a sparse DAG. Thus, the score function is 


F(W) = €(W;X) + AllW Ila. (5) 


To ensure acyclicity, they replace the constraint G(W) € 
DAGs in (1) with a equality constraint h(W) = 0, where W is 
a DAG if and only if 


h(W) = tr(e”°”) —d = 0, (6) 


where o is the Hadamard product, e” is the matrix 
exponential of a matrix M, and tr(-) is the trace of a squared 
matrix. The score function (1) then can be rewritten to a 
continuous optimization objective: 


ened IX —W*X||2 + Al|w]| 
w 2n p e 


s.t. tr(e”°”) —d=0. 


Equality-constrained optimization methods such as 
augmented Lagrangian method (ALM) [27] and asymptotic 
unconstrained method [20] can be used to optimize the 
objective. 


(7) 


In that we aim at study the sparsity problem, we simply 
use the nonlinear SEM proposed by [22] in the form of 


X = g,(W" gi (X) + Z), (8) 


with the ALM mentioned in [13], where g, and g, are 
nonlinear transformers. The nonlinear optimization objective 


is a combination of MSE and ?,-regularization with a DAG 
constraint: 


min (W, 0; X) +AllWIl,, 
st tr(eW™”) —d=0, 


where 0 is the parameter collection of g4 and gp. 


(9) 


III]. DEFECT OF THE ORDINARY SPARSITY CONSTRAINT 


We first theoretically investigate the phenomenon in Fig. 
1, then empirically analyze the problem by experiments. 


We start by discussing the problem of sparsity constraint 
with an interesting phenomenon in Fig. 1. There is a roughly 
light-colored vertical line on the 9" node of the estimated 
graph, indicating that all points except the dark-colored one 
on the 23" row are unexpected noises. Intuitively speaking, all 
noises would tend to be evenly distributed on the whole graph. 
However, that vertical line is rather counter-intuitive. 


Our theoretical hypothesis is as follows. Whenever the 
model makes wrong predictions, MSE in (9) traps into a local 
optimum that X takes all nodes on the column to fit itself, as 
Fig. 1 shows. Under the circumstances that MSE cannot 
escape the local optimum, 1.e., MSE is at local minimum, the 
optimization objective will try to minimize f4 -regularization 
and amplify the nonlinear transformers. Due to the form of 
sparsity constraint, it will minimize ?, -regularization by 
decreasing absolute maximum value of W rather than making 
the matrix sparse, because the noises on the column are 
effective for fitting X. As a result, the noise values of spurious 
edges will not be eliminated and take a relatively large ratio. 
We take a look at the optimization objective (9) first and then 
show some evidences. 


A. Theoretical Discussion about Optimization Objective 


We reorganize the unconstrained part regarding the j" 
node in optimization objective (9) to be 


n m 
1 
_ (k) (k) 
F(W) = md, X; — G2 3 Wii GG) (10) 
= J= 


alw; 


For brevity, we draw a comparison between two cases, 


perfectly fitted and misfitted, where the j node has no parents. 


Perfectly fitted case learns relation matrix well. May as 
l a a er 
well let the estimated matrix equals the true matrix, i.e., W, = 


W; = 0, then the j™ MSE part is 
1 
£;(W;,0;X) = 5 Var(Z;), (11) 
which means it only subjects to the variance of the j® node. 
For that the £4 -regularization term |W; I, = 0. So (10) equals 
1 
F(W) = 5 Var(Z;), (12) 


which only depends on random noise distributions. 


Assume the misfitted case reaches a local optimal solution. 
Let the estimated matrix not be consistent with the true matrix, 


L€; w # wr. If not well fitted, f; would take as much 


related information as possible into account, resulting the 
presence of nearly all nodes on the j™ column. Owning to the 
additive property of the optimization objective and 
nonlinearity, whenever MSE traps into a local optimum, if it 
cannot escape, the only thing F;(W) can do is to minimize £4- 
regularization |w; I, , and maximum components like 
nonlinear transformers like g, and g2. In this condition, the 
objective is 


min ty (Wp 0; X) + AIM, (13) 

where £;* is some constant local optimum with mutable W; 
and 0, thus (13) can be further written as 

min AW (14) 


Two strategies can be performed in minimizing WA ss 
decreasing the absolute maximum value in proportion or 
making the matrix sparse by discarding noise values. However, 
the ordinary sparsity constraint poses no restriction on the 
chosen of strategies, which may lead to absolute maximum 
value decrease of the matrix rather than making the matrix 
sparse, thus may not be able to reduce the ratio of the noises 


to the true values. 


B. Empirical Exploration of the Problem 


We use empirical evidences to support our theoretical 
hypothesis. Three experiments are conducted using the GAE 
model proposed by [22] and keeping default hyperparameters 
while setting number of variables to 50, i.e., the same 
experiment setting as is in Fig. 1. The first experiment aims at 
testing the hypothesis that when MSE reaches local optimum, 
the optimization objective will try to minimize the sparsity 
constraint. The second experiment tries to show that the 
objective tends to decrease the absolute maximum value of W 
as well as amplify the nonlinear transformers. Due to the 
problem of ordinary sparsity constraint, the noises cannot be 
reduced and take a relatively large ratio to true values in W, 
which is the last experiment illustrates. 


1) The first experiment: We compare the trends of MSE 
and the f,-regularization in Fig. 2. 

Fig. 2 denotes a typical case of the training curves of total 
loss, MSE and #,-regularization. We can clearly see that when 
total loss reaches the optimum at about 2500" epoch, MSE 
reaches its optimum at the same pace, while #4 -regularization 
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Fig. 2. Training curves of the total loss, MSE and f,-regularization. The 
red solid line is total loss, and the green dashed line and the blue dashdot line 
denote MSE and f,-regularization respectively. When total loss reaches 
optimum at nearly 2500" epoch, MSE reaches optimum at the same time, 
while f,-regularization keeps descending. 


keeps descending from 10° to 5 x 107?. It proves that when 
MSE reaches local optimum, the optimization objective will 
continue to minimize the sparsity constraint. 


2) The second experiment: We first show the absolute 
maximum weight of matrix W at different training epochs in 
Fig. 3, followed by a curve of the ratio between the median 
true weights and median false noises in training in Fig. 4, as 
well as the amplitude of the nonlinear transformers in Fig. 4. 
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Fig. 3. Visualization of the matrix weight at different epochs. We choose 
the 100", 500", 1000", 2000", 3000" and 5000" epoch for showing the 
changes of absolute maximum weight, which decreases rapidly. 
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Fig. 4. Training curve of the ratio between the median true weights Wp 
and median false noises Wp, as well as the amplitude of 0. The right vertical 
axis indicates how transformers amplify. After 2500 epochs, the ratio falls 
into stagnation, and the amplitude is larger than 1.0. 


Fig. 3 shows the absolute matrix weight in ordered epochs. 
It is clear that with the epoch increasing, the absolute 
maximum weight decreases rapidly. Fig. 4 is the training 
curve of the ratio between the median true weights and median 
false noises, as well as the amplitude of the nonlinear 
transformers. Shortly after the loss function reaches optimum 
at 2500" epoch, the ratio Wrp/Wep falls into stagnation, 
which indicates that although £} -regularization is still 
descending, optimization objective attempts to reduce the 
absolute maximum weight and amplify @ rather than make the 
matrix sparse by discarding noises. 


3) The last experiment: We show that the noises cannot 
be reduced and take a relatively large ratio in Fig. 5. 

We draw the minimum true weights and the maximum 
noise weights of each column at different epochs, for that good 
predictions rely on thresholding to distinguish true results and 
noises. We come to a conclusion from Fig. 5 that after the 
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Fig. 5. Visualization of the minimum true weights and maximum noise 
weights of each column in matrix at different epochs. We choose the same 
epochs as in Fig. 3 for showing the weight evolution. The noises take a 
relatively large ratio even at the 5000" epoch. 


2500 epochs, the noises are still unneglectable and take a 
relatively large ratio. Besides, the ordinary sparsity constraint 
is unable to converge quickly. 


We design three experiments to list evidences for 
supporting our theoretical hypothesis that when the 
component MSE cannot be reduced, the optimization 
objective may try to minimize the sparsity constraint by 
decreasing the absolute maximum value of W and amplifying 
the nonlinear transformers rather than making the matrix 
sparse, thus the noise values of spurious edges will not be 
eliminated and take a relatively large ratio, and will finally 
affect model performance. 


IV. A MODEL WITH PROPORTIONAL SPARSE PENALTY 


We present the model structure first, and then propose our 
proportional sparse penalty constraint. Training method is 
placed last. We named the model GAE-PSP. 


A. Graph Autoencoder Model with PSP 


The model uses graph autoencoder framework which can 
easily embed nonlinear SEM. We now connect the model in 
Fig. 6 with the equation in (8). 


Let g, and gz represent the 3-layer MLP encoder and the 
3-layer MLP decoder respectively, ignoring the random noise 
variable Z, the variables X, H and A can be written as 


H = g1(X), 
f = WTH, (15) 
X = g,(A). 


The corresponding optimization objective is to minimize 
the reconstruction error with the sparsity constraint. 
Concretely, the objective consists of MSE, the sparsity 
constraint and the acyclicity constraint. Thus we rewrite (9) as 


(16) 


s.t. tr(e”°”) —d=0, 


where 6, and 0, are the parameters of g4 and gz respectively, 
and Sy, is the proportional sparse penalty that discussed later. 
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Fig. 6. Visualization of the model structure. X is the n i.i.d. d-dimensional sample of m variables. W is the matrix representation of the causal structure of 
the m variables which is derived from the Hadamard product of a random initialized matrix and a masked off-diagonal matrix. Note the masked off-diagonal 
matrix is used for preventing variables predict them by themselves. Through a 3-layer MLP encoder, the model transforms X to the latent representation H, 
then computes H by matrix multiplication of H and W, and further gets X through a 3-layer MLP decoder. 


B. Proportional Sparse Penalty 


is not effective 


Ordinary sparsity constraint, 1.e. ZA p 
enough as we have discussed in Section III. A plausible reason 
for it is that 2, -regularization does not always make the matrix 
sparse. To solve the problem, we come up with the idea that 
the sparsity constraint should always focus on relative 
proportions rather than on absolute maximum values. We 


propose the constraint as 


Sw = lo (am) 


where W /max(|W]) is the proportional value of W, o and @ 
for controlling the sparsity are the sigmoid function and the 
hyperparameter respectively. With higher @, more small 
noises will be projected to the boundary of o, which means 
noises receive more penalty. 20(-) — 1 shifts the result in 
parentheses to the range between -1 and +1. 


l (17) 


1 


With Sy, the optimization objective computes the sparsity 
by the proportional weights without suffering from absolute 
maximum value decreasing, and imposes more penalty on 
small noises given higher y, which helps to make the matrix 
sparse. 


C. Training 


In order to solve the equality constraint optimization 
problem, we apply ALM which introduces a quadratic penalty 
to (16): 


1 Pe 
Ly (W, 01,62,@) = = |X — X]|_ + ASw iis 
+ah(W) +E AW), 


where h(W) = tr(e”°”) — d, æ is the Lagrange multiplier 
and p is the penalty parameter. We can solve the optimization 
problem by iteratively updating æ and p with Algorithm 1. 


V. EXPERIMENTS 


In this section, we provide a comprehensive evaluation on 
the proposed model, and investigate the performance of our 
proposed PSP constraint. We compare it against three recent 
gradient-based methods, DAG-GNN [24], NOTEARS-MLP 
and GAE [23]. We implement the proposed model based on 
PyTorch [28] and follow the default hyperparameters found in 
the three models except that we set sparsity constraint 
coefficient to 0.01/d* , œ in (17) to 2.0 and graph 
thresholding to 0.2. 


Algorithm 1 Train GAE-PSP 


Input: Observations X, initial values a) and Po, hyperparameters 
y and f, p’s max value Panax, and number of iterations K. 

Output: Weighted matrix W*. 

]: for number of iterations k + 1 E€ K do 

2: while p**1 < P x do 

3: Update W*+1, af**, 05+! by minimizing L,(W, 4,02, a) 
4: if h(W**1) > yh(W*) do 

5. pk*t e Bpkt 

6: else 

7: break 

8: end if 

9: end while 

10: qkt1 a ak + oXtth(w**2) 

11: end for 


We experiment on scalar-valued (d = 1) synthetic dataset 
using the same dataset setup as in [22]. Concretely, the ground 
truth DAGs are generated from the Erdds-Rényi model with 
four graph sizes m € {10,20,50,100} and expected node 
degree 3. We then assign random edge weights from a uniform 
distribution in (—2,—0.5) U (0.5,2). Sample X is generated 
3000 1.1.d. samples by sampling the nonlinear model 


X = 2sin(W! cos(X + 1) + 0.5- 1) 
+(W* cos(X +1) + 0.5-1) (19) 
+Z, 


which is the same as [23], under two types of random noise, 
Gaussian noise and Multivariate non-normal random noise by 
[29]. 


We choose Structural Hamming Distance (SHD), True 
Positive Rate (TPR) and G-Score which are widely used in 
this area. SHD measures the similarity between two graphs 
and smaller is better. TPR measures the ratio between the 
predicted true edges and all true edges and larger is better. G- 
Score is defined as 


max(0, (TP — FP)) 


C Score= 
a TP+FN. ” 


(20) 


and larger is better as well. 


A. Main Results 


TABLE I. , TABLE II. and TABLE III. list experimental 
results of the compared methods on SHD, TPR and G-Score 
respectively. For different number of variables m, we conduct 


10 experiments for each model with different seeds, and 
compute its mean and standard deviation. 


We can find that the performance of our proposed GAE- 
PSP consistently ranks best among other methods, and scales 
well under different number of variables. For SHD results, our 
model has lowest mean values among others, and the 
variances are much smaller than GAE. However, the ratio 
between variances and mean values in GAE-PSP is larger than 
that in DAG-GNN, it is probably because DAG-GNN fixes 
initial state of the relation matrix and introduces random noise 
representation which makes it more stable, and NOTEARS- 
MLP follows the principle of independent causal mechanisms 
(ICM) [30] which helps exclude spurious edges. For TPR 
results, GAE-PSP performs the best, which means it can 
accurately identify correct edges. GAE has the competitive 
performance as GAE-PSP. For G-Score, GAE-PSP 
outperforms other methods. In particular, there is a great 
performance improvement of GAE-PSP compared to GAE, 
which may be that GAE gets good TPR at the cost of 
predicting spurious edges, while GAE-PSP eliminates them. 


TABLE I. RESULTS ON SHD 
Gaussian Noise 
peer m=10 m=20 m=50 m=100 
DAG-GNN 12,144.72 24.5+6.62 64.349.43 | 135.7419.22 
NOTEARS-MLP* | 5.72.54 13.5+4.30 55.0+7.42 - 
GAE 6.6+13.09 | 20.0+56.84 | 26.0+153.38 | 31.3+65.12 
GAE-PSP (ours) 1.844.62 4.0410.55 3.3418.90 4.7+6.86 
Model Multivariate Non-normal Random Noise 
m=10 m=20 m=50 m=100 
DAG-GNN 13.4+5.19 25.946.95 | 68.1410.58 | 138.7420.72 
NOTEARS-MLP?*| 5.32.58 12.7+4.57 51.4+8.50 - 
GAE 10.8+15.96 | 26.9462.59 | 58.1+261.33 | 60.4+178.35 
GAE-PSP (ours) 2.344.14 6.0412.84 2.3+7.44 4.6+7.16 
* We generate 1000 samples for NOTEARS-MLP for all experiments 


because of large time consume in training. 
a. Experiment for NOTEARS-MLP on 100 variables takes too much time, 
thus, we fail to get the SHD, TPR and G-Score values. 


TABLE IL. RESULTS ON TPR 
Gaussian Noise 
MOESA m=10 m=20 m=50 m=100 
DAG-GNN 0.52+0.14 0.44+0.10 0.29+0.05 0.31+0.05 
NOTEARS-MLP | 0.65+0.10 0.58+0.09 0.37+0.70 - 
GAE 0.83+0.26 0.92+0.21 0.98+0.09 0.97+0.11 
GAE-PSP (ours) | 0.93+0.20 0.94+0.18 0.98+0.06 0.97+0.03 
Multivariate Non-normal Random Noise 
Model 
m=10 m=20 m=50 m=100 
DAG-GNN 0.55+0.12 0.43+0.10 0.28+0.05 0.31+0.05 
NOTEARS-MLP | 0.65+0.14 | 0.61 0.10 0.4+0.07 - 
GAE 0.79+0.29 0.90+0.19 0.97+0.11 0.98+0.08 
GAE-PSP (ours) | 0.91+0.20 0.91+0.25 0.98+0.06 0.97+0.04 
TABLE III. RESULTS ON G-SCORE 
Gaussian Noise 
More m=10 m=20 m=50 m=100 
DAG-GNN 0.19+0.20 0.18+0.12 0.16+0.09 0.14+0.08 
NOTEARS-MLP | 0.59+0.18 0.54+0.13 0.34+0.11 z 
GAE 0.71+0.39 0.82+0.36 0.87+0.26 0.83+0.30 
GAE-PSP (ours) | 0.89+0.27 0.89+0.26 0.9740.11 0.97+0.04 
Model Multivariate Non-normal Random Noise 
m=10 m=20 m=50 m=100 
DAG-GNN 0.15+0.20 0.13+0.12 0.12+0.08 0.12+0.09 
NOTEARS-MLP | 0.63+0.13 0.57+0.11 0.37+0.10 - 
GAE 0.57+0.43 0.74+0.41 0.85+0.30 0.80+0.33 
GAE-PSP (ours) | 0.83+0.29 0.83+0.32 0.97+0.10 0.97+0.05 


B. Comparative Experiments 


We conduct three experiments on our proposed GAE-PSP 
corresponding to those on GAE in subsection B of Section II 
to show superiority of our proposed PSP constraint. We use 
same experiment setting as in Fig. 1 to Fig. 5. 


— Loss 
=-=- MSE 


Metrics 


—-- ,-regularization 
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Epochs 


Fig. 7. Training curves of the total loss, MSE and f,-regularization. The 
red solid line is total loss, and the green dashed line and the blue dashdot line 
denote MSE and f,-regularization respectively. Total loss reaches optimum 
at nearly 1500" epoch, and MSE reaches optimum at the same time, different 
from that in Fig. 2, ,-regularization converges as well. 


1) The first experiment: We compare the trends of MSE 
and the f,-regularization in Fig. 7. 

We can see clearly in Fig. 7 that total loss, MSE and £4- 
regularization converges at the same time at the 1500" epoch. 
It proves that when MSE reaches local optimum, our 
optimization objective (16) has also minimized the sparsity 
constraint. 


2) The second experiment: We first show the absolute 
maximum weight of matrix W at different training epochs in 
Fig. 8, followed by a curve of the ratio between the median 
true weights and median false noises in training in Fig. 9. 

Fig. 8 demonstrates that after the objective converges, the 
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Fig. 8. Visualization of the matrix weight at different epochs. We choose 
the same epoch as in Fig. 3 for showing the changes of absolute maximum 
weight, which keeps almost unchanged after converges. 


absolute maximum weight of matrix W keeps stable as well. 
We can see that the absolute maximum weight stays 
unchanged in the last three subplots. Fig. 9 is the training 
curve of the ratio between the median true weights and median 
false noises. After the loss function reaches optimum at the 
1500" epoch, the ratio Wrp/Wpp keeps increasing, which 


— Wrp/Wrp 
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Epochs 


Fig. 9. Training curve of the ratio between the median true weights Wp 
and median false noises Wpp. The ratio keeps increasing, which implies that 
the objective keeps reducing noise values. 


implies that the optimization objective attempts to make the 
matrix sparse by discarding small noises, and that 1s what we 
expect the sparsity constraint should do. 


3) The last experiment: We show that the noises in our 
GAE-PSP model are small enough to be neglected in Fig. 10. 
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Fig. 10. Visualization of the minimum true weights and maximum noise 
weights of each column in matrix at different epochs. We choose the same 
epochs as in Fig. 5 for showing the weight evolution. The noises downgrade 
fast and are neglectable after nearly 1500 epochs. 


We come to a conclusion from Fig. 10 that the noises are 
neglectable after the optimization objective converges as the 
last three subplots show, and that our proposed objective (16) 
converges faster than (9) comparing the last four subplots 
between Fig. 10 and Fig. 5. An interesting phenomenon is that 
noises at the 9" column both in the estimated graph of Fig. 1 
and in the last subplots of Fig. 5 disappears in the last subplots 
of Fig. 10, which indicates that our proposed PSP constraint 
may be able to mitigate abnormal noises. 


In summary, we demonstrate that the PSP constraint can 
help eliminate spurious edges and largely improve model 
performance. The main results show that our proposed GAE- 
PSP outperforms other models both on the Gaussian and non- 
Gaussian datasets. And the comparative experiments show 
evidence that our proposed optimization objective (16) may 
work properly by reducing unexpected noises, thus making 
the matrix sparse. As a result, our proposed PSP constraint can 
solve the problem arising from ordinary sparsity constraint to 
a certain degree. 


VI. CONCLUSIONS AND FUTURE WORK 


In this paper, we observe the phenomenon that there are 
sometimes light-colored vertical lines in the estimated relation 
matrix. We investigate it and find it rises from the joint effect 
of the property of the optimization objective and the nonlinear 
transformers of the model. In particular, we theoretically and 
empirically find that when the model predicts spurious edges, 
i.e., it traps into a local optimum, the optimization objective is 
prone to minimize the absolute maximum value of the sparsity 
constraint and amplify the nonlinear transformers, rather than 
makes the matrix sparse. Inspired by that, we developed a 
proportional sparse penalty constraint named PSP, which aims 
to mitigate the phenomenon and improve model performance. 
We demonstrate its effectiveness by our proposed GAE-PSP 
model on the experiment comparing with other three models 
and verify its reliability on three comparative experiments. 
Future works include investigating the availability of PSP on 
other nonlinear SEM forms, studying the approach for better 
mitigating wrong predicted noises, and devising a way to 
reduce the sensitivity of graph thresholding. 


The code is at https://github.com/ThinkNaive/GAE-PSP. 
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